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Abstract 


Free electron laser theory is developed from the Maxwell and Lorentz force equations; 
the properties and characteristics of the laser are reviewed. The wave equation is solved 
for the fundamental Gaussian mode, and higher-order modes in Cartesian and cylindrical 
coordinate spaces, yielding expressions for the complete and orthogonal basis sets of Hermite- 
and Laguerre-Gaussian beams. Motivated by the evident inclusion of higher-order modes 
in free electron laser simulations, a tool is developed for the higher-order (in particular 
Laguerre-Gaussian) modal analysis of simulated free electron laser beams. 
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Chapter I. 


The Free Electron Laser 


In 1971 Stanford University physicist John Madey proposed the development of a coherent 
light sonrce based on the radiation emitted by relativistic free electrons: a free electron 
laserd In 1976, his team amplihed a CO 2 laser beam using relativistic free electrons from 
the Stanford Superconducting Linear Accelerator. They later observed coherent infrared 
radiation at a wavelength of 3 /rm in an oscillator conhguration. These achievements gave 
birth to the contemporary free electron laser described below. 

A. System 

The two main components of a free electron laser are an accelerator which imparts kinetic 
energy to a beam of free electrons, and an undulator also called a wiggler which extracts 
energy from the beam in the form of coherent radiation. There are a variety of FEL designs; 
all contain these two essential elements. Varying among individual conhgurations are the 
methods by which the electron beam is created, transported, then dumped, and how the 
radiation is amplihed. Some FEL’s are amplihers, some are oscillators. The oscillator 



Figure 1.1: Oscillator FEL 


conhguration is shown in Figure 1.1. Electromagnetic radiation is stored within a resonant 
cavity, and the increased optical power is outcoupled through a partially-transmitting dielec¬ 
tric. In the ampliher conhguration (Figure 1.2), spontaneously emitted light or light from a 
seed laser is amplihed in a single pass through the undulating electron beam. 

U. M. J. Madey, “Stimulated Emission of Brehmsstrahlung in a Periodic Magnetic Field”, Journal of 
Applied Physics, Vol. 42, No. 5, pp. 1906-1913, April 1971. 
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Electron Beam 



Figure 1.2: Amplifier FEL 


The source of free electrons is usually a cathode within an electron gun. Whether con¬ 
figured as an oscillator or amplifier, the electron beam itself may be either recirculated or 
dumped after exiting the undulator. Recirculation allows the recovery of beam energy back 
to the accelerator, enhancing overall system efficiency, with the added benefit of dumping 
less energetic electrons, producing less background radiation. In the RF recovery configu¬ 
ration, electrons return to the accelerator 180° out of phase with respect to those entering 
from the beam source. Figure 1.3 shows the entire system for a recirculating free electron 
laser oscillator. 


i 

recycled electrons 


electron beam 



optical 
cavity 
: mirrors 


n: 




beam 

dump 




source superconducting 

accelerator essential for 
recycling electron energy 


wiggler for conversion 
of electron energy 
into light 


Figure 1.3: Recirculating free electron laser oscillator. 


Attributes 

The free electron laser’s most striking difference from conventional lasers is also its chief 
advantage. The absence of a vulnerable medium such as a crystal or fluid chemical matrix 
means that an FEL can achieve extremely high intensity without damaging the laser itself. 
High wall-plug efficiency (~10%) is predicted for high power FELs with recirculating electron 
beams. High reliability is also a characteristic: currently operating FELs run continuously 
and reliably for weeks. 
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Free electron lasers are continuously tunable. The relativistic Lorentz factor 7 of the 
electron beam, the undulator period Aq, and undulator RMS magnetic held strength Rrms 
determine the laser wavelength A according to the resonance condition^ 


A 


Ao(l + i^') 

272 


where the undulator parameter K is dehned as 


K 


2711710 ^ 


( 1 . 1 ) 


( 1 . 2 ) 


Here e is the electron charge, m the mass of the electron, and c the speed of light in vacuum. 
The relation (1.1) suggests that the undulator wavelength and magnetic held as well as the 
energy of the electrons^ may be adjusted within an operating laser to yield light over a 
range of wavelengths. Free electron lasers have been designed to produce radiation from 
microwaves to X-rays, and they are continuously tunable over a smaller range within these 
regimes. In addition, the short-pulse nature of the radiation may enhance propagation of 
optical power through the atmosphere, making the laser particularly suited as a weapon. 


B. General Theory 

Free electron laser theory addresses the motion of the oscillating electrons, the propagation 
and dihraction of the resulting radiation (the optical mode), and the subsequent interaction 
of that radiation with the electron beam itself. We forego a discussion of diffraction until 
the next chapter, and proceed with a description of the electron motion, the optical mode, 
and gain. 

1. Electron Motion 

The trajectory of a relativistic charged particle in external helds is governed by the Lorentz 
force equation® 

^ (1.3) 

^At resonance in an FEL, one wavelength of light passes over one electron in one undulator period. 

^Tuning via electron beam energy is generally considered more difficult than via undulator geometry. 
Indeed, significant changes in electron energy require adjustments to guiding and focusing magnets along 
the entire length of the accelerator. 

^P. Sprangle, J. R. Pehano and B. Hafizi, “Propagation of Intense Short Laser Pulses in the Atmosphere”, 
Phys. Rev. E, 66, 046418, 2002. 

^P. Sprangle, J. R. Pehano, A. Ting, B. Hafizi, and D. F. Gordon, “Propagation of Short, High-Intensity 
Laser Pulses in Air”, J. Dir. Energy 1, 73, 2003. 

®We shall use T for proper time; r will be used for dimensionless time; cgs units are used throughout this 
chapter. 


3 



where = ( 7 c, 7 u) is the four-velocity of the particle and 


paf} 


/ 

0 

~Ex 

-Ey 

-E, 

\ 


Ex 

0 

-B, 

By 



Ey 

B. 

0 

-Bx 


V 

E, 

— By 

Bx 

0 

/ 


(1.4) 


is the field-strength tensor, whose components are given by the periodic magnetic held of 
the undulator and the electric and magnetic components of the radiation held. For the 
propagation of a plane wave along the 2 :-axis of a helical undulator,^ the undulator and 
radiation helds are taken to be 



= i?o(cos sin/co^:, 0 ), 

(1.5a) 

E,. 

= F'(cos'^, — sin'^, 0 ), 

(1.5b) 

Br 

= F'(sin'^, cos'^, 0 ), 

(1.5c) 


where Bq and ko = 2 ti/\q are the amplitude and wavenumber of the undulator held, E = 
E{x,y, z,t) is the optical held amplitude, and the phase is 'ip = kz — ut + (f){z,t), where 
k = ui/c is the optical wavenumber, and (l>{z,t) is the optical phase. The held-strength 
tensor is then 


pa(3 ^ 


/ 0 —E COS Ip Esinip 

E cos '^0 0 

—Esinp) 0 0 

V 0 — Bo sin fco^—-E cos'0 Bocosk^z^Esinij) 


0 \ 

Bo sin koz+E cos ip 
— Bq cos kQZ—E sin -0 
0 


Inserting (1.6) into (1.3) gives the four coupled diherential equations 
dU^ 


dT 

dU^ 

Hr 

dU^ 

hr 

dU^ 

Hr 


=- (U^Ecosp— U'^E sin pp) , 

me 

= -r U^E cos p) — (Bq sin koz + E cos p)] , 

me 

= —— \—U^E sin p; — {—Bocoskoz — F^sinV')] , 
me 

=- {~Bo sin kQZ — Ecos p) — {Bq cos kQZ + E sin p)] . 


Making the substitutions 


E = cp\ 

im dt 

dT = —, 
7 


( 1 . 6 ) 


(1.7a) 

(1.7b) 

(1.7c) 

(1.7d) 

( 1 . 8 ) 

(1.9) 


^In a helical undulator, consecutive magnetic fields are rotated 90° from each other causing the electron 
to follow a corkscrew path. The derivation for electron motion in a linear undulator is more complicated. 
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where 7 = (1 — /3 • /3) we can rewrite (1.7) as 

6 

— =- E cos'ip - f3y sin'ip) , (1.10a) 

at me 

^ilPx) = -—[EcospJ-(5z{Bosmkoz + Ecos^p)], (1.10b) 

at me 

d 6 

— - [—E simp — Pz {—Bq cos koz — E simp)] , (1.10c) 

CLl TTXC- 

d 6 

-piilPz) = - [E {P^ cos ^p - Py simp) + Bo{Px sin koz - Py cos koz)] . (l.lOd) 

(XT' T1X(^ 

The magnetic field does no work, and indeed (1.10a) states that the electron’s energy is 
affected only by the coupling of the radiation electric field to the electron’s transverse motion. 
It is convenient to rewrite this equation as 

7 = -—/3-E. (1.11) 

me 

By defining a transverse velocity /3j_ = + Pyj = {P^, Py, 0), we may combine (1.10b) and 

(1.10c) into a single expression 


d(^ l3j_') 6 

- - -=- [E (1 — Pz) (cosip, — sin'^, 0) + PzBq (— sin k^z, cos koz, 1 

at me 

which for highly relativistic electrons {Pz —> 1) can be approximated by 


d(7/3j 

dt 


~- [PzBq (— sin koz, cos k^z, 0)1. 

me 


Equation (1.13) is readily integrated: 

J d( 7 / 3 ±) = 

7/3± = 

P± = 


eBf 


0 


me 

eBp 

me^ 


eBn 


kp'ymc^ 

eBpXp 


Pz{— sin kpz, cos kpz, 0)dt, 
^(— sin kpz, cos kpz, 0)pX, 
cos kpz, sin kpz, 0) + C 
2 (cos kpz, sin kpz, 0) + C. 


2'K'yme^ 

We eliminate the constant of integration by assuming perfect injection, so that 

- -w 

Using (1.2) and the fact that for a helical undulator i?rms = Bp, we can write 

K 

P3± = -(cos/co:^, sin/co;^, 0) 

7 


( 1 . 12 ) 

(1.13) 

(1.14) 

(1.15) 

(1.16) 

(1.17) 

(1.18) 
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or 


jSx = -cos koz, 


(1.19) 

( 1 . 20 ) 


Substituting (1.18) into (1.11), we have 

eKE 

7 = -(cos/co^:, sin0) ■ (cos'^, — siu'^, 0) 

■ymc 

= - -(cos/co^: cos'^ — sin/co^: sin'^) 

■ymc 

= - cosikoz + ip). 

ymc 

Meanwhile, substituting (1.19) and (1.20) into (l.lOd) gives us 


-y-iy/3z) =- E[ - cos koZ cos ^jJ-\ - sinkozsimjj 


f K K 

+Bo -cos koz sin koz H-sin koz cos k^z 

V 7 7 

cK 

yPz + yPz = - [E (cos koz cos ifj — sin koz sin i/j) 

ymc 

— Bq (cos koZ sin k^z — sin koZ cos kQz )], 

cKE 

y[3z = - cos{kQZ + Ip)- yp^. 

ymc 


Substituting (1.21) into (1.24), we have 


cKE eKE 

y/3z = - cos(koz + ip) - cos(koz + ip)/3z, 

ymc ymc 

eKE 

y(3^ = (l-pz) - cos{koz + ij), 

ymc 

eKE 

/3^ = (1-/3^)— — cos{koz + ^p). 
y^mc 


( 1 . 21 ) 


( 1 . 22 ) 


(1.23) 

(1.24) 


(1.25) 

(1.26) 


(1.27) 


It is useful at this point to collect the electron terms in koz + '^ by identifying the electron 
phase 

C = {k + ko)z — ut, (1-28) 

referenced to the combined optical and undulator helds. This leads to 


koz + Ip = C -3- (p- 


(1.29) 


c = (k+koty-uj 

at 

= {k + ko)/3zC- u. 


(1.30) 
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Solving for (3^ gives 


^ = 

$z = 


c 

{k + ko)c 

C 


{k + ko)c' 

The expression for (3^ rnay be simplified as follows: for laser light we require 

2tic 


u = 


A 




Meanwhile, the change in electron phase over the course of the undulator® is 


: 27rc 9 

C ~ —— ~ 10 s 
Ao 


In other words, since typically Aq 3> A, then cu 3> C and (1.31) becomes 

B _. 

(/c + ko)c 

Using (1.29), (1.32), and (1.35), we may now rewrite (1.27) as 

C eKE 


(1.31) 

(1.32) 

(1.33) 

(1.34) 

(1.35) 


U 


COs(C + 0), 


{k + /co)c 7^mc [ {k + ko)c\ 
eKE 

C = -7—; [{k + ko)c - u] cos(C + 0) 


'y^mc 
eKEko 


cos(C+ •)>). 


(1.36) 


(1.37) 


We can normalize the evolution of the electron motion to the length of the undulator by 

o 

dehning a dimensionless time r = ctjL. We adopt the notation C, = d(/dT and we write 


c = 


eKEkpK 

2TTeKEK 


COs(C + 0) 
COs(C + 0) 


2TiNeKEL 

-5-^-COs(C + 0), 


(1.38) 


where N = L/Xp is the number of undulator periods. Dehning the dimensionless optical 
held 

I I _ 2nNeKEL 

|a| =-^-r—, (1.39) 


®To achieve gain, the electrons must “bunch” in phase as they traverse the undulator. This is discussed 
further in the next two sections. 
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we have the pendulum equation 


C= |a| cos(C + <t)) 


(1.40) 


Figures 1.4 and 1.5 show the results of a simple simulation of an operating FEL. The electron 
trajectories in phase space are indeed of a pendular nature: the paths are elliptical (closed 
paths, as in Figure 1.4 - the pendulum swings back and forth), and there exists a separatrix, 
a boundary outside of which change in electron phase ( does not reverse sign (open paths, 
as in Figure 1.5 - the pendulum goes “over the top”). These hgures will be discussed more 
in the following sections. 



Figure 1.4: Electrons moving in phase space, within the separatrix {i/q = 0). 
Although there is some evolution of the optical phase, there is no significant 
gain on resonance. (The shading of electrons darkens with increasing t.) 



Figure 1.5: Electrons moving outside the separatrix = 2.6). There is gain 
due to favorable bunching of the electrons in phase space. 






2 . Evolution of the Optical Mode 


The fields generated by the oscillating electrons must satisfy the inhomogeneous Maxwell 
equations 





c 


( 1 . 41 ) 


where 




/ 0 —E cos Ip Esm-pj 

E cos '^0 0 

—Esinp) 0 0 

\ 0 —E cos Ip Esin Ip 


0 \ 

E cosp! 
—E sin Ip 

0 / 


( 1 . 42 ) 


is the radiation part of (1.6), and J" = (cp, J) is the electron beam 4-current. We seek a 
4-vector potential A" = (<h, A) satisfying 


^/3 ^ 


( 1 . 43 ) 


and the Lorenz condition (9 q,A" = 0, for which (1.41) becomes the wave equation 

((9"A^ - 

^ ^ c 

A 'jr 

dad'^Af^ 

c 

c 

where O = (9^(9" is the d’Alembertian operator. A suitable choice® is 


(1.44) 


which gives 


A“ = 


/ 0 \ 


^ • / 

— SlU'^ 

k 

— cosip 
k 


\ 0 / 




d E . , d E , ^ 

— — sm V -f- — — cos V = 0 
ox k oy k 


( 1 . 45 ) 


( 1 . 46 ) 


®This choice is suitable because it corresponds to a helically polarized wave. A general expression for A“ 
in terms of the electron 4-velocity is available in the form of the Lienard-Wiechert potentials (see pp. 661-665, 
J. D. Jackson, Classical Electrodynamics, Wiley, New Jersey, 1999). Our choice is made with the benefit of 
hindsight and is not unique. 
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and, using u = kc, 


<924° 

dA^ 




= 0, 

dxo 

dxo 

dA^ 

dA^ 

E d 

dxo 

dxi 

kcdt 

dA^ 

MO 

E d 

dxo 

dX2 

kcdt 


(1.47a) 

(1.47b) 

(1.47c) 


satisfying (1.42). We now write (1.44) as two second-order partial differential equations 

1 E . , E . , Att ^ , 


1 E , E , 47r ^ 


For E{z, t) and 0 ( 2 :, t) slowly varying, that is 


d‘^E , dE 
^ ^Ih ’ 


d‘^E dE d‘^(j) d(j) d‘^(j) dcj) 


we keep only the boxed terms in 

1 d^E . , 2 dEd4> 


2 dEdcp , 2 dE 

-— COS lb - cos w 

ujc dt^ c dt 


E 3“^ (t> , E /d(P 

-— cos w -— 

LOC dt^ LOC V 3t 


2E dcp ioE 

-sini/>- smip— 

c dt c 


d‘^E , dE , 2c dEd<p 

-—^sin'0 -2^—cos'ip -—^cosf/; 

' az^ az uj oz^ 


cEd'^E , luE , ^90 , ioE /d4>\‘^ 

-:^cos'0H-sin'0+ 2E — simp H- I — I simp 

UJ dz^ c dz c \dzJ 


(1.48b) 


(1.49) 


(1.50a) 


leaving 


1 d'^E , 2 dEdtp , 2dE , E d'^p , E f dpY 

LOC ot^ UJC at^ c at UJC ot^ ujc \dt J 

2E d<t> cd'^E dE 2cdEd<t> 

+- COS'Ip - cosip - —cos Ip2 - simp H-^sin-^ 

c dt c UJ dz^ dz UJ dz^ 


cE d^E 


sin xp H-cos-i/^H- 2E—-cosip + 

dz^ c dz 


ujE / d(f) 


(1.50b) 


d(j)\ , fldE dE\ Ae ^ . . 

fldd) (90\ fldE dE\ 47r / .n 

Dehning J_i_ = {J^, Jy,0), this can be written 

fldE dE\ , , , , fldd) dd)\ , , , ^ 47r , , 












and identifying the orthogonal unit polarization vectors 

€i = (—cos'^, sin'^, 0) 

€2 = (sin'^, cos'^, 0) 


we can write 


The operator 


dE dE\ ^ ^ 

^ + = 2n},-e„ 




d d 


d dz d 


dt ^ dz dt dt dz 


(1.53a) 

(1.53b) 


(1.54a) 


(1.54b) 


(1.55) 


is the rate of change in coordinates traveling with the evolving light wave, and we write 


dE 2nL 

dr c 

d(j) _ 271L 

t,— — - J±-^2- 

dr c 

The transverse current is the sum of individual electron currents, that is 


(1.56a) 

(1.56b) 


J± = -ec^/ 3 i ( 5 (x - fj), 


(1.57) 


or, using (1.18), 


(5(x — rj) (cos koz, sin k^z, 0). 


(1.58) 


Equations (1.56) are then 


dE 277 KeL 


(5(x — rj)(— cos k^z cos ip + sin kozsinip)^ (1.59a) 


M 27TKeL^^^ 7/ , . , . , 

E— = - > o(x — rj(cosA;o 2 :sin'^ + sm/co^cos'^), 

rin- ^ f ^ 


(1.59b) 


which may be written 


277 KeLup 


dr 7 Jo 

^dcp 277KeLne 

dr 7 Jo 


COs(C + 0) —, 
sin(C + 0)^, 

277 


(1.60a) 


(1.60b) 
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where Ue is the electron number density, and individual electrons are identihed by their 
initial coordinates in phase space (Co, ^o)- All electrons are assumed to have the same initial 
energy z/q- These may be combined into 


dE d(j) 


2TiKeLne 

7 

2n KeLrie 
7 



i sin(C + (f>) 


dCo 

27r’ 


Multiplying both sides by gives 


dr 7 Jo 2n' 


A compact form of this expression is found by identifying 


(1.61) 

(1.62) 


(1.63) 


a 

j 


a e 


Sn^K^e^L^Nue 


(1.64) 

(1.65) 


where |a| is given by (1.39). Letting an overbar indicate averaging over all initial electron 
phases, we have 

( 1 . 66 ) 



It is clear that there is no change in optical held {da/dr = 0) for j = 0. For gain, the 
integral further requires that the initial phase distribution of electrons average to less than 
zero. Populations of electrons for which 7^ 0 are said to be bunched in phase, and 
the manipulation of electron phase distributions to achieve favorable bunching is a critical 
concept in FEL design. 


3. Interaction and Gain 

We now have expressions for the evolution of the electrons in the undulator (1.40) and 
the light they create (1.66). To show how this light is amplihed, we introduce a solution to 
(1.40) for weak helds (|a| tt) and low current (j < tt). It is an expansion in oq, namely^® 

C('r) = Co + t'oT - ^[cos(Co + vqt ) - cos(Co) + z/orsin(Co)] H -, (1-67) 

^0 

o 

where z/ = C, Oo is the initial optical held, and the initial phase is taken to be zero. Thus^^ 

^{ t ) = 1^0 + — [sin(Co + i ^ qt ) - sin(Co)] 

^0 

^2 r 

+ -|l - t[cos(2Co + 2z/or) - cos2Co] (1.68) 

^0 L ^ 

+ cos z/qt - 1 - z/qt sin Co cos(Co + 1 H- 


B. Colson, “Free Electron Laser Theory”, doctoral dissertation, Stanford University, 1977. 
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We can use (1.67) and (1.68) to plot the phase space portrait of populations of electrons over 
their travel time through the undulator, but a more accurate picture is achieved numerically 
using a computer simulation, such as the one used to generate Figures 1.4 and 1.5. 

The amplihcation of the optical mode is proportional to the energy loss of the electron 
beam. Leaving (1.68) behind for the moment, we can examine the change in individual 
electron phase velocity v by differentiating (1.28): 


dr 


d 

^ = -^[{k + ko)z - ut] 


L 

c 


/ 7 1 \dz 

(fc + fco)— -uj 

at 


= L[{k + ko)f3z - k], 


(1.69) 


which, assuming /c 3> /co, may be approximated 


u = Lk{P,-l), (1.70) 

that is 

Au = LkAf3z. (1.71) 

This is not a surprising expression: change in electron phase velocity is certainly propor¬ 
tional to change in its longitudinal velocity, as the Lorentz force felt by the particle as it 
passes over each magnet in the undulator depends on this velocity. Moreover, proximity to 
the resonance condition—most certainly a function of longitudinal electron velocity!—will 
ultimately determine the gain. Near resonance, we may write (see Equation 1.1): 

^ Ao(l + iF2)’ 

thus 

= MiTlo) 

1 + K^ 


(1.72) 


(1.73) 


To put this expression in terms of electron energy, we must re-examine the relationship 
between Pz and the relativistic 7 . In our helical undulator, where (1.18) implies 


the dehnition 




2 



1 


(1.74) 


(1.75) 
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becomes 


1 


7 = 


'1 - /3i - (3^ 


(1.76) 


which may be written 


/ 

7 


1 + 7^2 


(1.77) 


l-(3i = — 
7 


and so, taking the differential of each side, 

13 = 


A7(1 + 772) 


(1.78) 


(1.79) 


A/3,= 


A7(1 + 772) 


Plugging this into (1.73) we have 


Az/ = 


47rAA7 


(1.80) 


(1.81) 


At this point we may return to our approximation —>■ 1 for relativistic electrons, and we 

write 


Az/ = 47rA—. 


(1.82) 


For an electron of energy Ee = 'jmc^, the energy change as a function of the change in 
electron phase velocity is then 


AEe = ^A.. 

47rA 


(1.83) 


Now, the number of electrons in a volume VJ, of electron beam is UeVb. The change in electron 
beam energy AEb corresponding to 14 is thus 


, ^ rie'ymc ^ ^ 


(1.84) 


where we have averaged the change in phase velocity for individual electrons. Meanwhile, 
the optical mode energy E^ in a volume 14 is simply 


E^ 

Eo = —14, 

47r 

where E is the optical field strength. We now write the optical field gain as 


(1.85) 


( 1 . 86 ) 


Uejmc^Vb^ 

Uejmc^E-^ 


(1.87) 
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where the hlling factor F specifies that fraction of the optical mode coincident with the 
electron beamd^ Using our dehnitions for the optical amplitude (1.39) and beam current 
(1.65), this can be simplihed to 

G = ( 1 . 88 ) 

To hud an expression for Az/, we return to our weak held expansion of z/(r) (1.68) and 
integrate over all possible initial electron phases 

Az/ = z/ — z/q 


r*27r 




( Oo 


r*27r 


sin(Co + vqt ) - sin(Co) 


27r 


+ ^ ^ [ cos(2Co + 2z/or) - cos 2Co] 


+ cos z/qT - 1 - z/qT sin Co cos(Co + f H- 


WCo- (1.89) 


Many of the terms average to zero, and so, to second-order in oq, we have 

2 cos z/qT — 2 + z/qT sin(z/or) 


Az/ ^ Og 


2z/3 


(1.90) 


Plugging this into (1.88) and dividing through by jF, we have 


G 

2 cos z/gT — 2 + z/gT sin(z/gr) 


z/3 

‘^0 


(1.91) 


which describes the PEL gain along the undulator. The hnal FEL gain at r = 1 (the end of 
the undulator) is shown in Figure 1.6. The salient feature of this hgure is its antisymmetry, 
in particular its behavior close to z/g = 0,^^ which indicates that energy may be either 
transferred from the electron beam to the optical mode (0 < z/g) or vice-versa (z/g < 0). 
Earlier (see p. 7), we remarked that electrons must bunch in phase space to transfer energy 
to the optical mode. For a functioning, useful FEL, (1.91) narrows that requirement to a 

^^For electron beam and optical mode volumes both spanning a longitudinal length I, the filling factor is 
defined as a ratio of cross sectional areas: 

Vb _ T^rlX _ A 

Vo irrlyF Ao 

^^It is interesting to examine why (1.91) does not blow up at z/g = 0- After all, the presence of z/q in the 
denominator is certainly suggestive of singular behavior at the origin. The answer is found by expanding 
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Figure 1.6: The FEL weak field gain curve plotted at the end of the undulator 
(r = 1). Energy is either gained by the optical mode or lost, depending on 
initial phase velocity t'o- 


population of electrons possessing an average phase velocity just to the positive side of the 
G/jF axis in Figure 1.6. 

To illustrate, we refer once again to Figures 1.4 and 1.5. These show the results of a 
simple simulation of FEL gain for an electron beam with z/q = 0 and z/q > 0. In the hrst case 
(Figure 1.4), nearly the same amount of energy gained by the optical mode from electrons 
with decreasing phase velocity is simultaneously lost to electrons whose phase velocity is 
increasing. In the second case (Figure 1.5), where z/q > 0, electrons bunch in the desired 
region of phase space: they are predominantly losing energy to the optical mode. 


the numerator about vq = 0, where 

2cosr'o'r — 2 + r'oTsin(r'oT) 


2 - 


., 4.^4 


2 2 , ‘'0 ‘ 




12 


-2 + iyoT { vqt - 




12 


In the vicinity of the origin, then, at r = 0 (1-91) becomes 

G 1 

— = —r'n, 
jF 12 


confirmed by the linear nature of the t'o ~ 0 region of Figure 1.6. 
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The introduction of the gain equation completes our preliminary discussion of free elec¬ 
tron lasers. For convenience, we review the important conclusions of this chapter: the two 
principle equations of basic FEL theory, and an expression (derivable from the hrst) for gain 
in weak helds. 


Pendulum Equation 

00 

c= \a\ cos(C + 0) 

Wave Equation 
a = —j 

Gain Equation 

(weak fields) 

G 2 cos z/qT — 2 + uqt sin(z/o'r) 

jF ^0 


(1.40) 


( 1 . 66 ) 


(1.91) 
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Chapter II. 


Fundamental Mode Optics 


The free electron laser is unique among high power lasers in that it operates primarily 
in the fundamental mode. In contrast, high-power solid state and chemical lasers deliver a 
beam comprised of many higher-order mode (HOM) components. Over a given propagation 
range, these ROMs lead to more diffraction loss than, for example, a purely Gaussian beam. 
While such a beam characterizes the ideal laser, we will show in the next chapter that 
beams obtained from simulations of FELs appear to contain a few ROMs, and it is our goal 
to construct a diagnostic tool to evaluate the mode content of a given laser beam whose 
amplitude and phase are known at every point. 

Such comprehensive knowledge of beam makeup is possible in the laboratory with much 
effort, but a simulation provides all such information a priori and allows a comprehensive 
understanding of the system, and the ability to alter that system, without the high cost 
of a laboratory experiment. These types of simulations have been benchmarked against 
actual experiment and have accurately predicted experimental results. Modal analysis, on 
the other hand, has not been benchmarked, and this analysis provides a starting point for 
such a comparison. 

A. Propagation - Analytic Description 

Although we intend to approach the FEL numerically, we are interested in an analytic 
solution to the wave equation in order to understand the nature of an evolving wavefront 
inside a laser cavity. For the moment, we consider the optical mode independently of its 
interaction with the electron beam. 

1. The Parabolic Wave Equation 

The homogeneous part of (1.44) can be written 

( 2 , 1 ) 

Expressions solving this equation describe waves propagating in free space (or other homo¬ 
geneous, sourceless media). Using the same vector potential as in the previous chapter, we 
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have 

= (gin^, cos'^, 0), (2.2) 

k 

where e = (—i, 1, 0) is the polarization vector, a = kz — ut is the carrier wave phase, and 
£{x,t) = i?(x, is comprised of an amplitude E and phase 0 now with an explicit 

transverse dependence (remember that before we used E = E{z,t) and 0 = (f){z,t)] see 
page 10). To begin substituting this into the wave equation (2.1), we take two derivatives in 
2 ; and t: 


A(x, t) = 3? 


T(x, 


dA 

dz 

d^A 

dz"^ 

dA 

d^A 

'W 




3ft <! ^ fSf - 2iu^ - u^S 1 e*“e 

dt 


k \ dE 


(2.3) 

(2.4) 

(2.5) 

( 2 . 6 ) 


As before, for slowly-varying E{'x,t) and 0(x, t), variations in amplitude and phase arising 
from diffraction in the direction of propagation are small compared to those of the wave 
itself, that is 


d‘^E 


, dE 

dz‘^ 

•C 

k — 

dz 


d'^E 


dE 

dE 

< 

Co’—— 

dt 




dz"^ 




k^ 

^dz 


dd'cj) 

dE 





(2.7) 


These negligible terms comprise the second longitudinal space and time derivatives of 8, 
which we eliminate to obtain 


d^S d^8 


dx"^ dy"^ 


2ik— - k^8 -^ I —2iuj— - uES 


dz 


dt 


e*“e = 0. 


( 2 . 8 ) 


This simplihes to 


V^8 + 2ik[^^ + ~ 
dz c dt 


e*“e = 0 , 


where we have used u = kc and introduced the transverse Laplacian 

<92 92 


= 


+ 


dx'^ dy"^ 


Multiplying by ue the scalar equation is 

. , [98 dS , 

cV^£ + 2 ,k (— + I - 0, 

and we may once again use (1.55) to write this as 


1„9 ik d . ^ 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


( 2 . 12 ) 
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With propagation now expressed in terms of r, we would like to make the diffraction term 
dimensionless also. The 1/e radius wq of the initial, fundamental mode of the laser is given 


by 



Zq\ 

TT 


(2.13) 


where zq is the Rayleigh length, or the characteristic distance over which the beam cross 
section increases. Dehning the dimensionless transverse coordinates as 



so that 


_ k <92 

dx"^ 2L dx^ ’ dy"^ 


and adopting the dimensionless transverse Laplacian 


k 92 


92 92 


we now write the dimensionless parabolic^ wave equation 



(2.14) 

(2.15) 

(2.16) 

(2.17) 


2 . Solving the Dimensionless Parabolic Wave Equation 

To solve (2.17), we attempt a trial solution of the Gaussian form2 

^{x,y,T) = £oPir)exp 


( x2 + y2 \ 

V <i{r) )’ 


(2.18) 


where p(r) and g(r) are as yet undetermined complex functions describing the evolution of 
amplitude and width over dimensionless time. We further impose the initial condition 


£(x,y,0) = To exp 




(2.19) 


^The terms paraxial and parabolic may both be used to describe (2.17). It is paraxial in that it accurately 
describes wave behavior close to the axis of propagation. It is parabolic in that it fits into a class of partial 
differential equations of the general form 


^ d'^u ^du du 

^ ^^dxdy ^ ^ ^ 


+ Fu = R{x,y), 


where 

= AC - = 0. 

This class of equations includes the Schrodinger equation, with which the parabolic wave equation shares 
many properties. In contrast, the original wave equation is hyperbolic, that is AC — < 0. 

^A spherical wave is a known solution to the wave equation. A Gaussian wavefront may be described as 
a paraxial spherical wave, an approximate solution to the wave equation. It is also, in fact, an exact solution 
to the paraxial wave equation. 


A B 
B C 
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where Zq = Wq is the dimensionless Rayleigh length and the central held amplitude is Sq = 
£^(0,0,0). Substituting (2.18) into (2.17), we take the appropriate derivatives in x, y, and r 
(letting and suppressing the argument r): 


d^S 

dx"^ 

d^S 

dy'^ 



-) Sope-^"/\ 

<lj 

Sope-^-"/\ 


( 2 . 20 ) 

( 2 . 21 ) 


d£ (Idp r‘^dq\ _^- 2 / 


( 2 . 22 ) 


and write 

(^--) 

V r qJ 


E^pe-^^'^ + 







0. (2.23) 


which may be written^ 


For this equation to hold for all r (that is, 
must be identically zero, thus 

dq 

dr 

dp 

dr 

It is easy enough to solve (2.26) for g(r): 


(2.24) 

(2.25) 

the expressions in parentheses 

(2.26) 

(2.27) 

go + IT- (2.28) 


Dividing through by and collecting terms, we have 

1 . fldp P dq' 


q^ q 


“ T + M r“7r + 


pdr q^ dr 


p dr qj 
for all X and y). 


= i. 


We are not surprised to hud that it is complex, as it must describe not only the evolution of 
the transverse extent of the wave, but also the evolution of the phase, which must oscillate. 
From our initial conditions, we hud 

go = zo, (2.29) 

also anticipated. Using this result and substituting (2.28) into (2.27), we have. 


dp ip 

dr zo + ir 

P 

izo - r’ 

^To get from (2.24) to (2.25), multiply through by —iq^, then collect terms in and q^. 


(2.30) 
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which we solve for p(r): 


Iq izq - r 


= \n{izo - t) 


= hiiizo) — ln(i5o ~ '^) 


= In 


r ’ 

l + i- / 
zo/ 


(2.31) 


p(t) =-5-. 

l + i — 
zo 


(2.32) 


What does it mean for p(r) to be complex? Our solution to the parabolic wave equation is 




=^ exp - ' . (2.33) 

l + i— V Zq + it J 
zq 

Once again letting we elucidate this expression by making the following 

dehnition: 

3 = 1 + ip, (2.34) 


where 


and we write 


£{x,y,T) = —exp - — 

3 V zoi 


So ( 3 

— exp -^ 

3 V zo 33^ 


= -exp -( 1 -*^) 

3 L ^o33 


y2 ~2 

■V-h ip- - 

2 ^ 033 * 2 ( 033 =' 


^0 f \ (■ 

— exp -T—^ exp ip -—^ 

3 V zoU* J \ ^o33 


Now, the complex number 3 has modulus 


(2.34) 


(2.35) 


(2.36) 


I 3 I — \/l + P^, 


(2.37) 
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and argument 
so we may write 

and (2.36) becomes 


arg 3 = arctan /i, 

g —2 arctan/i 


3 1/1 + ’ 


Sn 


: exp 


Defining the dimensionless beam width 




^—2 arctan (i 


W 


= 'J 5o33* = V5o(l +/i2) = , zo(l + ^], 


we arrive at 


S{x,y,T) = f - ~o/ x ) exp () exp ( -i arctan ^ ) , 

w(r) V M;2(r)y \4 + ^ J \ ^ 0 / 


where 


Wo = w( 0 ) = v^- 

This is the solution for the dimensionless parabolic wave equation. Using 


0(r,r) = 


r^r 
+ r 2 


T 

— arctan —, 
zo 


we write finally 


— C ^0 ^-r^/w'^{T)^i4>(f,T) 


S{x,y,T)=So—e 


(2.38) 

(2.39) 

(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 

(2.45) 


B. Properties oe the Gaussian Beam 

We have identified a physical basis for each of the factors in (2.45); these are reviewed 
in Table 2.1. Figure 2.1 illustrates a beam diffracting from a waist at r = 0 in the ±r 
directions, identifying the important parameters used in the Gaussian beam expression. We 
can study the behavior of the beam described by this expression by looking at its variables 
in their limits. As we have required, at r = 0, 

S{x,y,0)=Soe-^"/^\ (2.46) 

As r increases, so does w: the beam expands. For a given r, the overall value on the right 
hand side of (2.45) decreases with increasing r: the field strength attenuates. 

Figure 2.2 shows a Gaussian beam propagating through free space. The expansion and 
attenuation of the beam are evident. Also evident is that the Gaussian beam retains its func¬ 
tional form: there is a scaling of the overall structure, but the wavefront remains Gaussian 
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Table 2.1: Symbols Used in the Dimensionless Parabolic Wave Equation 


Symbol 

Definition 

Description (dimensionless) 

X 


orthogonal transverse coordinate 

y 


orthogonal transverse coordinate 

r 

1 71 r 

transverse coordinate 


ct 

time 

T 

L 

w 

h (' ^ 

beam width 

Wo 

w(0) = \/io 

beam width at beam waist 

Zo 

zo 

L 

Rayleigh length 


after propagation. This behavior is unique to unmixed, individual component solutions of 
the paraxial wave equation, but as we will show in Chapter IV., the equation allows solutions 
beyond just the fundamental Gaussian. Moreover, we will show that an arbitrary laser (for 
example, a laser with a top-hat or annular or just very “messy” cross-sectional wavefront 
structure) experiences the inclusion of higher-order modes—that is, higher-order solutions 
to the wave equation—as a loss mechanism^ during propagation. 

The above claim demands to be addressed, but to further motivate our discussion of 
higher-order modes, we first present the results of past research in the field of high-power 
free electron lasers. 


^This is not loss in the sense of absorption or scattering, but rather with respect to the amount of energy 
intended to reach a target area. 
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Figure 2.1: A diagram of the optical mode showing the relevant (dimension¬ 
less) physical meanings of the symbols in Table 2.1. 



Figure 2.2: The beam described in Figure 2.1 on the interval 0 < r < 1. 
The initial Gaussian wavefront (left) propagates through a uniform sourceless 
medium (center), diffracting into a wider Gaussian of lower peak amplitude 
(right). 
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Chapter III. 


Motivation: Free Electron Laser 
Simulations 


A computer simulation of laser beam propagation through a sourceless medium was used 
to generate Figure 2.2. Such a simulation implements a relatively simple numerical method 
for solving the parabolic wave equation. More relevant simulations (to wit: the inclusion of 
the radiation source—the electron beam!) are the work of the Free Electron Laser Research 
Group. Naval Postgraduate School FEL simulations have been in use for hfteen years and 
have been useful tools in the construction of new, more powerful FELs. Recently, ever 
increasing computing speed and power and advances in simulation technology (better codes) 
have resulted in greater hdelity, and more robust analysis of the evolving FEL system. The 
inclusion of expanding coordinates and methods of characterizing optical beam quality are 
among the improvements found in the following work, completed over the last two years. 

A. Multiple Modes 

Under optimum conditions, the free electron laser operates in the fundamental mode. In 
order to create the laser beam, the gain medium—a beam of free electrons—passes through 
the undulator, and thus some portion of the optical cavity. The complicated boundary 
conditions introduced by the presence of a region of swiftly moving charges within the re¬ 
gion of propagation argue strongly for the use of numerical methods. But one might guess 
that the simple evolution of an initially Gaussian wavefront like that in Eigure 2.2 is made 
more complicated by the presence of the electrons, and in fact simulations that include the 
electron beam show that an initially Gaussian FEL beam wavefront can assume an other- 
than-Gaussian transverse structure. 

1. Multi-Mode Simulations of a Short-Rayleigh Length FEL 

Simulations offer the ability to test variations in free electron laser design parameters 
such as undulator length, electron bunch charge, and Rayleigh length, before such designs 
are realized. For a weapon-class FEL, two main objectives are maximizing output power and 
minimizing the possibility of damage to resonator optics. In 2003, the FEL Group simulated 
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the performance of a proposed 1 pm high-power FEL with a Rayleigh length an order of 
magnitude shorter than typical extant operating FELsd The three-dimensional simulation 
(in X, p, and r) used self-consistent FEL pendulum (1.40) and wave (1.66) equations to 
follow the evolution of weak optical helds to steady state, allowing for the development of 
transverse structure in the wavefront. For simulations of this type, the value of interest is 
the extraction rj, dehned as 

extracted optical power 
initial electron beam power 


More specihcally, there are two expressions for rj in basic FEL theory:^ one for the low gain 
regime (j < 1), 

( 3 - 2 ) 


and one for high gain [j 3> 1), 


V 


2iV’ 

SN ' 


(3.3) 


Both are approximate, as is the dehnition of the cut-off between high and low gain. 

In the simulations, design parameters are independently varied to determine their opti¬ 
mum values with 77 as a gauge. The simulations must determine the optimum resonance 
parameter 

iyo = L[{k + ko)/3^o-k], (3.4) 


where L is the undulator length, k is the optical wavenumber, ko is the undulator wavenum¬ 
ber, and (3zo is the initial electron velocity (see Chapter 1 for details of these parameters). 
An actual FEL operates at the optimum z/q, Ending it by automatically varying A = 2T[/k 
until a fixed value of z/q dominates by mode competition. 


2 . Results 

Each individual simulation run results in output like that shown in Eigure 3.1. A series 
of runs are inspected to determine the optimum z/q for a given parameter, say charge q per 
electron bunch. The plot of extraction rj versus q at each optimum z/q shown in Figure 3.2 is 
thus the result of more than one hundred simulation runs. Figure 3.2 demonstrates the steady 
increase in extraction with increasing bunch charge. The simulation assumes a constant 
longitudinal emittance 

El = Ay/fe, (3.5) 

where Ay is the rms spread in the electron energies and If, is the rms electron bunch length. 
In a real FEL (and in later simulations), emittance increases with bunch charge, reducing 
extraction as the beam quality degrades. But the effect is estimated as small here.^ 

^W. B. Colson, et al, “Multi-mode Simulations of a Stiort-Rayleigh Length FEL”, Proceedings, 25*^ 
International FEL Conference, Tsukuba, Japan, September 2003. 

^W. B. Colson in W. B. Colson, C. Pellegrini, A. Renieri (Eds.), Laser Handbook, Vol. 6, North Holland, 
Amsterdam, 1990. 

^For a concise description of experimental FEL performance as it relates to particular phenomena, such 
as electron beam emittance, see Chapter 3 of B. W. Williams, Jefferson Lab Free Electron Laser 10 kW 
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Figure 3.1: Output from a 2003 simulation. Each output figure displays 8 
panels showing, clockwise from top left: (1) a plot of the mode evolution over 
many passes through the undulator, (2) a cross section of the laser beam, (3) 
a list of input and calculated parameters, (4) a plot of output beam power 
evolution, (5) a plot of gain evolution, (6) the final phase space portrait of the 
electrons, (7) an historical plot of electron phase space evolution, and (8) a plot 
of the mode evolution over the final pass through the undulator. 


The results of the 2003 runs, which were presented at the International Free Electron Laser 
Conference in Tsukuba, Japan, included determinations of optimum number of undulator 
periods {N = 14, see Figure 3.3), optimum normalized beam waist radius {a = 0.12, see 
Figure 3.4), and constraints on normalized Rayleigh length (Jq < 0.06, versus the typical 0.5 
in use at the time, see Figure 3.6) and mirror transmission losses (>50%) for the proposed 
laser. 

Figure 3.3 reflects two important FEL design concepts. The first is that too few undulator 
periods result in a non-functioning FEL, indicated by the rapid drop off in extraction for 
N < 10. This is evident from our dehnition of the dimensionless current density 


j 


Sn‘^K‘^e^L‘^Nne 

STT^K^e^XlN^Ue 


(3.6) 


Since gain G is proportional to current j, we ultimately have 


G (X j (X N^. 


(3.7) 


Upgrade - Lessons Learned, Technical Report NPS-PH-05-001, Naval Postgraduate School, Monterey, May 
2005. 
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Figure 3.2: Plots of PEL performance as a function of individual design 
parameters, in this case bunch charge q, are each the result of many simulation 
runs. 


Meanwhile, the dependence of extraction on number of undulator periods initially obeys 
(3.2), that is 

(3.8) 

countering the trend due to increasing gain. In the high gain regime, extraction appears to 
have no dependence on N: 


V 


8N 


m 






(3.9) 


1287^mc^ ’ 

and indeed the curve levels off once maximum gain is achieved, around iV = 15 for this laser. 

Figure 3.4 illustrates the effects of two important FEL concepts: filling factor F and 
(normalized) transverse emittance given by 


£n = inOb, 


(3.10) 


where 7 is the electron beam energy, Vb is the rms beam radius, and 6b is its rms angular 
spread. In general, low emittance is desirable, and the value can only increase as the beam 
travels around the FEL transport lattice. Considerable pains are taken to suppress emittance 
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Figure 3.3: Plot of extraction 77 as a function of number of undulator periods 
N. While the curve is relatively flat beyond N = 15, a guiding principal in 
PEL design is using just enough undulator periods to get a working, efficient 
laser. 


growth, but it is never zero, as evidenced when the electron beam is focused at the center of 
the undulator. Since £„ must at least be conserved, the narrower the waist of the electron 
beam, the greater the angular spread. Too narrow a waist, and some electrons fall outside 
the optical mode at the ends of the undulator due to a large angular spread, as in the top of 
Figure 3.5. On the other hand, the electrons in a wider beam interact with a lower average 
optical held in the middle of the undulator, as in the bottom of the hgure. Increasing F by 
increasing Ah may result in some electrons falling outside the optical mode, as illustrated in 
Figure 3.5. These electron still wiggle, but the radiation they produce is not reinforced by 
the optical resonator. 

The important result from Figure 3.6 is that for a small increase in extraction achieved 
by increasing the Rayleigh length zq, there is a considerable increase in optical intensity on 
the resonator mirrors. This is of particular concern at present. Since vulnerability of the 
mirror limits the allowable power, material limitations make the value of mirror intensity a 
limiting factor in the development of weapon-class free electron lasers. 
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Figure 3.4: Plot of extraction rj and optimum initial phase velocity t'o as 
a function of normalized electron beam waist radius a for fixed transverse 
emittance e„. There is an optimum overlap between the electron beam and the 
optical mode. 



Y 





Figure 3.5: The overlap of the electron beam and the optical mode. A nar¬ 
rower electron beam means larger emittance e„ and a lower number of inter¬ 
acting electrons, yet a wider beam creates an interaction region with a lower 
average optical field. 
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Zg, Normalized Rayleigh Length 

Figure 3.6: Plot of extraction 77 and mirror intensity as a function of nor¬ 
malized Rayleigh length Zq. For longer Rayleigh lengths, a small increase in 
extraction is accompanied by a large increase in (potentially harmful) optical 
intensity on the resonator mirrors. 
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B. Optical Mode Distortion 


The simulation used in the 2003 runs was adequate to the task of determining the results 
described in the previous section; scientists at Thomas Jefferson National Accelerator Fa¬ 
cility (JLAB), where the proposed laser was being developed, made use of the information 
provided, and the consideration of shorter Rayleigh length conhgurations was taken up as a 
serious endeavor."^ 

For the following year’s runs, a number of improvements were made to the simulation code. 
A faster fast Fourier transform (FFT) algorithm® was employed resulting in a 10-fold speed 
improvement. A more accurate propagation method involving next-nearest neighbors was 
utilized, and an expanding coordinate system that follows the rapidly diffracting wavefront 
made for greatly improved efficiency. 

1 . Optical Mode Distortion in a Short Rayleigh Length FEL 

A proposed high power FEL has the parameters shown in Table 3.1. As before, the 
simulation has the ability to follow multiple and arbitrary transverse optical modes, as they 
interact with the electron beam and bounce back and forth in the optical cavity, including 
mirror transmission and edge losses.® Output from the improved simulation is shown in 

Table 3.1: Parameters for a Proposed High-Power FEL 


Parameter 

Symbol 

Value 

electron beam energy 

Eb 

80 MeV 

bunch charge 

Q 

400 pC 

bunch length 

k 

0.15 mm 

peak current 

i 

800 A 

normalized emittance 


3 mm-mrad 

beam radius 

Tb 

60 pm 

undulator periods 

N 

22 

undulator wavelength 

-^0 

2.36 cm 

undulator length 

L 

52 cm 

peak magnetic held 

B 

0.7 T 

undulator parameter 

K 

1 

optical cavity length 

S 

18 m 

Rayleigh length 

zo 

6 cm 

resonator quality factor 

Qn 

4 

optical wavelength 

A 

1 pm 


^At this writing, an experiment is in place at JLAB involving a resonator configuration utilizing a Rayleigh 
length on the order of z = 0.2. 

^M. Frigo and S. G. Johnson, “FFTW: An Adaptive Software Architecture for the FFT”, Proceedings, 
1998 International Conference on Acoustics and Signal Processing, p. 1381, 1998. 

® J. Blau, et al, “Optical Mode Distortion in a Short Rayleigh Length Free Electron Laser”, Proceedings, 
International Free Electron Laser Conference, Trieste, Italy, September, 2004. 


34 



Figure 3.7. The displayed data (top center of the figure) clearly exhibits an other-than- 
Gaussian wavefront for the particular parameters of that simulation run. 


rn, mwEinowrr cvouition 



Figure 3.7: Output from a 2004 simulation. New visualizations of the optical 
mode wavefront accompanied improved speed and diagnostics. 
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2 . Results 


When the wavefront pictures provided by the simulation were included in the plots of the 
parameters being varied to hud optimum extraction (Figures 3.8 and 3.9), it became clear 
that under certain configurations,'^ a FEL operates in higher-order modes. It is interest- 


5 


4 





3 

2 

1 



Figure 3.8: Plot of extraction g as a, function of dimensionless Rayleigh length 
zq and wavefront pictures generated by the simulation. 



Dimensionless Mirror Separation, 


Figure 3.9: Plot of extraction g as a function of dimensionless mirror separa¬ 
tion Tniir and wavefront pictures generated by the simulation. 


ing to note that while varying the Rayleigh length zq (Figure 3.8) and mirror separation 
Tmir (Figure 3.9) evoke a relatively flat response from the FEL in terms of extraction, the 
new boundary conditions imposed by these changes appear to lead to a change in modal 
composition of the beam. 

^For example, under certain Rayleigh lengths or distances between resonator mirrors. Remember that 
we said the FEL automatically finds the optimum vq, so each data point in Figures 3.8 and 3.9 represents 
an operating FEL, with associated pictures indicating the shape of its optical wavefront. 
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C. The Need eor Higher-Order Mode Analysis 


The radially symmetric, generally papilloid shape of the higher-order mode wavefronts 
observed in the 2003 and 2004 simulation results is suggestive of the geometry of the FEL 
resonator system and the boundary conditions imposed by an optical wave propagating along 
and around a coaxial column of current (see Figure 3.10, below). Under the assumption that 



Figure 3.10: Wavefront suggesting coincident reinforcement of the optical 
mode by the electron beam. This is the general shape we see in many FEL 
simulations, which display the beam structure at the high reflector end of the 
resonator. The appearance of such structure demands a detailed analysis of 
the modal content of the outcoupled beam. 


the fundamental Gaussian mode is the ideal mode both for lasing efficiency and atmospheric 
propagation, a simple index of beam quality® would be sufficient. The observation from 
simulations that higher-order modes occur within the resonator demands a more robust 
characterization of beam quality. In other words, we wish to know more precisely which 
modes comprise our beam. 


®One such measure of beam quality frequently encountered in the literature is given by 



where d is the observed beam diameter and 



is the diameter of a theoretical diffraction-limited Gaussian beam arising from the same initial conditions. 
For more on see Safnett, and SiEGMAN 1993. 
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Chapter IV. 


Higher-Order Modes 


An approach to the resolution of a given laser beam into its component modes is a process 
analogous to a Fourier decomposition of a function. In the same way that an arbitrary 
function^ may be represented by a series of sines and/or cosines, an arbitrary laser beam 
may be represented by a sum of orthogonal modes. Identifying the appropriate complete, 
orthogonal basis set is the challenge, and is the topic of this chapter. 


A. Hermite-Gaussian Beams 


The trial solution we used to solve the parabolic wave equation in §1. can be further 
generalized to allow for higher-order modes. We begin by noting that (2.45) can be expressed 
as a product of identically dependent functions of each transverse coordinate, that is 

S{x,y,T) = 

W 

— Iw'^ Q-y'^ 1^'^ ^iy'^Tl{zl+T^) ^-i&vcta.n{Tlzo) 

W 



arctan(r/io)/2g—gix^r/(ig+T^) 


X 



arctan (r/zo)/2g-y2/u>2 


gW(^o+'^^) 


= S^{x,T)Sy{y,T). 


(4.1) 


As such, we can simplify the process of hnding higher-order solutions by examining only one 
transverse dimension at a time. Writing the wave equation as 


4 dx"^ dr ) 




0 , 


(4.2) 


^The function must, of course, meet the Dirichlet conditions: on any finite interval, it must possess one 
minimum, one maximum, and no discontinuities. 
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and modifying our trial solution to give the complex amplitude function p an explicit trans 
verse dependence, 


we begin taking derivatives: 

dx^ 


d‘^p A^dp 


+ -^x 


d£ 


dx"^ q dx \q 


p 


SqC 


-x^lq 


TF + P 
ox q^ dr 


to give 


-F + 'IF _ \ = 0. 

4 dx' q dx q' 2q dr q' dr^ 


Dividing through by —jEqc and collecting terms in spatial derivatives of p gives 


iFp 4 ^ 

OF q dx 


1 + 1— -^x 


dr J q 


p + 4# = 0. 

or 


It is the similarity of (4.7) to the well known Hermite differential equation 

d^Hm dHm 

— 2u —^-h 2mHm = 0 


(4.3) 

(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 


dn? du 

that, after some algebra, leads to the following expression for the Hermite-Gaussian beam:^ 


£{xf) = 

w 


^-X^lzo gi(/>moG,T) 


w 


or in three dimensions 


Wo 


S{x,y,r) = So—Hm 

w 


'V2: 


X 


w 




'V2y' 


w 


^-r^lzo^ilpmn{f,T) 


(4.10) 


where 

~2 

T T T 

4>mn{r, t) = —r - {m + u + 1) arctau —, (4.11) 

Zo + r2 ^0 

and Hm and Hn are rrd^- and n*’^-order solutions to (4.8). These so-called Hermite-Gaussian 
solutions to the dimensionless parabolic wave equation constitute a complete set of orthog¬ 
onal modes with which one may construct an arbitrary laser beam wavefront. 

The properties of Hermite polynomials are well documented.^ The hrst ten polynomials 
are listed in Table 4.1. Some examples of the effect of the functions Hm and Hn on the beam 
wavefront is represented in Figure 4.1. These “burn patterns” illustrate the intensity £‘^ of 
the beam in the transverse plane, but display nothing about the phase! To fully characterize 
a beam, one must know both its amplitude and phase at every point, a challenge in a real 
experiment, but a matter of course in a simulation. 

^The (dimensional) functional forms of £i(i,y,T) as well as a derivation of the higher-order solution are 
available in A. E. Siegman, Lasers, University Science Books, Sausalito, 1986. We include the dimensionless 
forms of the r-dependent expressions for completeness. 

^See, for example, pp. 367-9, V. Barcilon, “Special Functions” in C. E. Pearson, Ed., Handbook of Applied 
Mathematics, Van Nostrand Reinhold, New York, 1983. 
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Table 4.1: Hermite Polynomials - Solutions to 


+ 2mHm = 0 




Ho{u) 

H,{u) 

H2{u) 

Hsiu) 

H,{u) 

H,{u) 

H,{u) 

Hjiu) 

Hs{u) 

Hgiu) 


1 

2 u 

4^2-2 
Su ^ — 12m 
16m^ - 48m2 + 12 
32m® — 160m® + 120m 
64m® - 480m^ + 720m2 - 120 
128m^ - 1344m® + 3360m® - 1680m 
256m® - 3584m® + 13440m^ - 13440m2 + 1680 
512m9 - 9216m^ + 48384m® - 80640m® + 30240m 



Figure 4.1: Some Hermite-Gaussian modes. The transverse electromagnetic 
mode (TEM) subscripts correspond to the values of m and n in the polynomials 
Hjyi and Hn • 


B. Laguerre-Gaussian Beams 


In cylindrical coordinates, the dimensionless transverse Laplacian is 


1^ / 

f df \ df) r2 ’ 


(4.12) 


and (2.17) becomes 

Ids f dS\ 1 d^S dS 

r — + -^ H- 

r dr \ dr J QQ^ dr 

We might guess that this coordinate system is more appropriate to the resonator geometry, 
whose eigenmodes are certain to reflect the circular shape of the mirrors. Following the same 
procedure as in the previous section, a generalized trial solution to (4.12) yields a differential 
equation similar to Laguerre’s differential equation 



M- 


d^Lp 

du^ 


+ (m + 1 - r) + pLp = 0. 
du 


(4.14) 
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Like the Hermite polynomials, the polynomials that solve (4.14) also constitute a complete 
set of orthogonal functions, the hrst ten of which are given in Table 4.2. 


Table 4.2: Laguerre Polynomials - Solutions to u 


(pLp 

dv? 


+ (1 


\ dLp 


H- pLp — 0 


Lo{u) 

Li{u) 

L2{u) 

L^iu) 

La{u) 

L^iu) 

L&{u) 

Lj{u) 

L^iu) 

Lq{u) 


1 

I — u 


1 — 2u + 4^2 
1 — 3u + — 4^3 

1 — 4m + + i^u'^ 

1 - 5m + 5m2 - 1^3 + _ ^^5 

1 _ 6m + f M^ - f M^ + |m^ - + ^^6 

i-7u + fu^-fu^ + fy-iu^ + ^^u^-^y 

1 - 8m + 14m 2 - f m 3 + ffM^ - ^m3 + ^M^^ - ^M^ + ^M® 

1 - 9m + 18m2 - 14m3 + ^M^ - |m3 + + ^^8 _ ^^9 


The Laguerre-Gaussian beam, given by 


S{r,e,T) 


Wo 
^Q — 
W 






(4.15) 


where 

~2 

0!”(r, 6. t) = J ^ r, — {2p -I- m -I- 1) arctan — -|- mO, (4.16) 

is a dimensionless solution to the parabolic wave equation in cylindrical coordinates.^ Beam 
propagation is shown for the first ten Laguerre-Gaussian modes in Figure 4.2. Evidently, the 
wavefronts retain their shape. Examining (4.15), we see there is no change in the number of 
zeros as r increases: there is no “crossing” of local maxima. 

^The subscripts m and n are standard in the literature for the Hermite polynomials. Similarly, m and p 
are often used for the Laguerre indices. At this point, we essentially leave behind our discussion of Hermite- 
Gaussian beams, so from now on m will denote the angular mode of the Laguerre-Gaussian solution. 
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Figure 4.2: The first ten Laguerre-Gaussian modes. Notice that Lq is simply 
a Gaussian beam. Notice also that all modes retain their general shape—in par¬ 
ticular the number of zeros—as they propagate through a sourceless medium. 


The orthogonality of the generalized Laguerre polynomials is weighted by “ 
normalization^ 

r u^e-^L^{u)L^{u)du = + ^ + 

Jo Q- 


with the 

(4.17) 


where 


r(a) = 


e ^du. 


(4.18) 


Jo 

We shall only consider those Laguerre modes with no azimuthal dependence, that is m = 0. 
For integers (which q + 1 surely is) 


r(n) = (n — 1)!, 


(4.19) 


and so, for our purposes, we shall write the orthogonality condition as 

Lp{u)Lq{u)du = 5pq. (4.20) 

Now, hnding the Laguerre-Gaussian (LG) “spectrum” of an arbitrary initial wavefront 
E{r, 0) with azimuthal symmetry requires us to calculate the coefficients Cp in the series 

E{r, 0) = ^ CpSp{r, 0) (4.21) 

p 

^See pp. 533-4 in D. Zwillinger, CRC Standard Mathematical Tables and Formulae, Ghapman & Hall, 
New York, 2003. 
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where 


^pix 1 0) ^oLp ^ I e 




zo 


-X/zo 


(4.22) 


solves (4.13) at r = 0. Notice that the factor in (4.15) becomes unity since 

X 0 0 

0°(r, 0, 0) = - (2p + 0 + 1) arctan-h 0 x 0 

^ z^ + 0 ^ ^ ^ zo 

= 0 . 

Expressing our arbitrary wavefront as a function of 


0 = 


2r^ 


and noting that® 




_ 2r^ _ 2r^ 

Qo - 2 ’ 

Zo 


we multiply both sides of (4.21) by an LG mode Sq{go) and integrate over qq: 

poo poo _ 

/ E{go)Sq{go)dgo = / ^CpSp{go)Sq{go)dgo 

Jo Jo p 

poo 

= Sq'^Cp e~‘^^Lp{go)Lq{go)dgo 
P 

= Sq CpSpq 


(4.23) 


(4.24) 


(4.25) 


- £^c 
‘-"0 + ’ 


®For example, using this new coordinate the Gaussian beam becomes 


(4.26) 
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where we have used (4.20). An arbitrary coefficient is then^ 


^2 I ^{Qo)^q{QO^dQQ 


■'0 Jo 

4 


^0^0 


‘2t‘^ 

~g-r2/zo^^ / __ I 
^0 


(4.27) 


Equation (4.27) is limited to r = 0. Even as we assert the usefulness of a modal analysis 
of the initial wavefront, we wish to analyze our beam at any point along its propagation, 
that is for any r. Beginning once again with (4.21), this time we make no assumptions about 
the value of r and thus w{t) and the r-dependent phase factor 


E{q,t) = ^CpSp{Q,T) 
p 


(4.28) 


Implicit in this generalization is the possibility of a complex Cp(r). Here is why. Suppressing 
the arguments g and r, let us multiply both sides of (4.28) by arbitrary mode S* and 


^We can normalize the quantity to unity by setting E{f) to the fundamental and calculating 


where we have used 


Co = 


Jq 

4 


^OZq 

4 






'0 


_ ^ £o 

£oZo 4 

1 


r[(n+ l)/p] 

pa{ri+i)/p ' 


with n = 1, a = 2/zo, and p = 2 (see p. 451, Zwillinger). Having thus identified 

CO 

and returning to our standard LG index p, a normalized Laguerre-Gaussian spectral coefficient is then 

'2f^' 


/^O , 


ZO 


E{f)df, 


with the carat indicating normalization. While this is a useful result, we ultimately intend to approach the 
calculation of LG-spectrum coefficients numerically, and a self-consistent normalization will be included in 
the computer program. The above is included for completeness. 
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integrate: 


poo poo _ 

/ E£*dg = / '^CpSpSgdg 

Jo Jo 


^o'y\cp / ^Lp{g)Lq{g)dg 

p Jo 

poo 

Sl^Cpe^iip-<i)^^^^^'^(rrzo) / e-<^Lp{g)Lq{g)dg 

p Jo 

c2 2i(p-ij)arctan(T/io)x 

Oq / ^ L-pC Upq 


- £^C 

^0 

and so, switching indices from q to our standard p, 


Cp = / E£*dg 


<72 


or explicitly 


~ /•lxj 

C„ = ^e*(2p+darctan(r/.-o) / Lp(g)E(g, T)dg 

JO 


^gi(2p+l)arctan(r/.-o) / r — E(h, r)dh, 

tow Jo V wV 


where once again 


3 = 1 + *^- 
2^0 


(4.29) 


(4.30) 


(4.31) 


(4.32) 


We will attempt to solve (4.31) numerically, for two reasons: (1) while we may know the 
amplitude and phase of a simulated FEL beam at every (hnite) transverse point, we will not 
generally have an analytic expression for the entire wavefront; and (2) even if we did know 
the functional form of E(r,r), the integral in (4.31) is difficult to perform for all but the 
most ideal beams. 

For a numerical analysis, we must hrst get the integral into a form the computer can 
understand. Letting the subscripts r and i indicate the real and imaginary parts, we rewrite 
(4.30) as 

1 /■“ 

Cp = j2 


C2 / (-^r “1“ djQ 




C2 / ^^p,rEi tSp^iEj' “h Sp_^iEi^ dQ 




C2 / (,^p,rEr Sp^iEj^ dQ “h ^2 / i^^p^rEi Sp^iEj-^ dQ 




^iJo 


(4.33) 



We may now split Cp into its real and imaginary parts, each of which may be integrated 
separately: 



1 

poo 

/ {^p,rEr + £p,iEi) 

Jo 

1 dg, 

(4.34a) 

III 

1 

^0 ■ 

poo 

/ (S E- - E ■E'' 
Jo 

1 dg. 

(4.34b) 


The contribution to the LG spectrum of an arbitrary beam by an individual LG mode shall 
be given by the modulus of that mode’s coefficient, 


Cp = \cp\ = + 

and will be normalized to the sum of all basis set coefficients: 


(4.36) 


C. Numerical Propagation of Higher-Order Modes 

The diagnostic tool to be included in free electron laser simulations numerically integrates 
(4.34) and reports Cp for each of the hrst ten LG modes via (4.35). Figure 4.3 is a sample 
output from bm.c, a self-contained code that resolves an input laser beam wavefront into 
component LG modes and reports the normalized LG coefficients as a bar chart beneath 
the wavefront picture both before and after propagation. In this case, the integration has 
correctly identihed the initial wavefront as an exclusively LG-2 beam: in the bar chart below 
the initial wavefront intensity plot, only the bar in the third position is present and goes 
to the top of the chart. To within a few percent error, the program has also determined 
the hnal wavefront to be LG-2, as we expect. The text in the middle of the output reports 
parameters used for that simulation as well as a “fractional energy on target” for a circular 
target of hxed radius. This is discussed in the next section. The basis sets used to analyze 
the input beam both before and after propagation are evaluated analytically and displayed 
at the bottom of the output for reference. 

The diagnostic is able to pick out multiple modes, as in Figure 4.4. In this case, we see for 
the hrst time the “crossing” of zeros in the propagation of the laser beam. Remember that 
“pure” LG modes do not exhibit this behavior: a given Laguerre-Gaussian mode possesses 
the same number of zeros throughout its propagation (see again Figure 4.2). Single modes 
also always exhibit a global maximum at r = 0, not always observed in mixed modes. 
This distinction between single- and mixed-mode higher-order beams immediately suggests 
a useful application for our tool, and it leads us to consider one of the difficulties in the 
experimental analysis of laser beam modes. 

It is not possible to determine the modal composition of a beam (function) if only its 
intensity (square) is known. As a simple illustration, let us consider the curves in Figure 4.5. 
It is not possible to determine whether the curve on the right was arrived at by squaring the 



(4.35) 
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Time Steps (ts) = 20000 (13888 required) 
Wndow Size (W) = 18.0G00G0 
"V^ist LcKation (tauw) = 0,000000 
Rayleigh Length (zO) = 0,500000 
Nunter of Bins (n) =300 
Target Radius (wni) = 0.450000 
Fractional Energy on Target = 0.107658 



C=0,003472 
eiTor=0,030366 


Figure 4.3: Output from the diagnostic tool bm.c. The basis set before and 
after propagation (evaluated analytically) is included at the bottom of the 
display for reference. 


top left or bottom left function. An attempt to determine the modal composition of a laser 
beam based on intensity information alone encounters the same problem. It is precisely the 
phase information that is lost in the process of squaring. 

Consequently, burn patterns can be misleading. Figure 4.6 shows the propagation of the 
same laser beam as Figure 4.4, but takes the cross sectional amplitude (or burn pattern) of 
the hnal wavefront earlier than before. Visually inspecting the initial and hnal wavefront 
amplitude plots (enlarged in Figure 4.7), one might be hard-pressed to conclude (or even 
guess!) that they were generated by the same laser beam. Yet the values of the LG coeffi¬ 
cients, as reported in the bar charts beneath the initial and hnal wavefronts in Figure 4.6, 
assure us that this is the case. Multiple modes interfere because each mode prescribes a 
unique time dependence via the factor 


^ gxp 


r^r 


r 




L ^0 


— i{2p + 1) arctan — 
^0 


(4.37) 


in (4.15). This evolving phase-front results in an evolving amplitude-front as modes cyclically 
cancel and reinforce each other, as evident in the propagation plot. 

The diagnostic tool performs well for even more complicated mixtures of higher-order 
Laguerre-Gaussian modes. In Figure 4.8, bm.c has correctly identihed the 3:5:2 mixture of 
the LG-0, -3, and -6 modes comprising the simulated beam. Figure 4.9 shows the analysis 
of a beam comprised of many LG modes. Once again, in both hgures we observe the zero 
crossings typical of mixed-mode beams, and the difference in the burn patterns before and 
after propagation. 
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Figure 4.4: A beam composed of two distinct Laguerre-Gaussian modes. The 
slight error in the bar chart pertaining to the LG coefficients for the final 
wavefront is a result of numerical integration error. 



Figure 4.5: When squared, the two functions sin x (top left) and sin(a; + tt) 
(bottom left) lead to an identical squared function (right). 


D. Power on Target 

Knowing which modes comprise our laser beam allows us to examine the relationship 
between physical properties of the beam and its composition. An example of this is the 
important consideration of power on target: how much of the beam’s initial energy is focused 
onto a target at the propagation distance. 

Choosing a circular target of an arbitrary radius such that almost all (93.5%) of the 
outcoupled energy is delivered to its surface when the beam is in the fundamental, Figure 4.10 
plots the energy delivered to the same target by beams in unmixed, higher-order Laguerre- 
Gaussian modes. 

Evidently, by themselves higher-order modes are not good choices for delivering power to 
a target at the center of the beam. As p increases, more and more energy is relegated to the 
“shoulders” of the central held (see Figure 4.11). Once in this outer region, light cannot cross 
back into the center without a change in the number of local held maxima—not observed in 
pure LG beams. 

Certainly, real lasers are unlikely to be comprised solely of single higher-order Laguerre- 
Gaussian modes, or even predictable mixtures. Figure 4.12 shows three theoretical lasers: 
one a Gaussian, one an even mixture of the Gaussian and LG-1 modes, and one a beam with 
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Figure 4.6: Each LG mode carries a different phase factor, causing periodic 
interference when multiple modes are combined. 



Figure 4.7: A combination of LG modes before (left) and after (right) prop¬ 
agation over a finite distance (enlarged from Figure 4.6). While the two wave- 
fronts are not necessarily identifiable by eye as manifestations of the same laser 
beam, the diagnostic assures us this is so. 


a “flat-top” cross section. The percentage of energy received by the centrally located target 
is shown. 

What the propagation of mixed modes such as the flat-top beam teaches us is that the 
inclusion of higher-order in the laser does not necessarily spell doom for the propagation of 
power to the target area. Certainly the trend in Figure 4.10 is not towards larger central held 
amplitude for higher values of p. But mixed modes, as Figures 4.6, 4.8, and 4.9 demonstrate, 
have the potential to get some of that power from the shoulders back into the center of the 
beam. In fact, we may conclude that any non-fundamental beam delivering more power 
to our test target than the Li-Gaussian is by definition comprised of multiple higher-order 
modes. Particularly for use as a weapon, if a laser beam is not in the fundamental, it may 
beneht from the mixing of higher-order modes. 
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Figure 4.8: A beam comprised of three higher-order modes. The recipe is 
three parts LG-0 (the fundamental Gaussian), five parts LG-3, and two parts 
LG-6. Once again, the burn patterns before and after propagation are difficult 
to reconcile visually. 
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Figure 4.9: A beam comprised of many higher-order modes. 



p 


Figure 4.10: EO and El are the energy contained in the beam and the en¬ 
ergy delivered to the target area, respectively. For unmixed Laguerre-Gaussian 
modes, laser power on target drops off significantly for p > 0. 
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Energy in the “shoulder” falls outside the target area. 



Figure 4.11: A target just large enough to receive nearly 100% of the energy 
in the fundamental will generally not receive energy from pure higher-order 
modes beyond the first zero-point. 
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Figure 4.12: The propagation of (a) a Gaussian beam; (&) a beam comprised 
of equal parts Lq- and Li-Gaussian modes; (c) an initially “flat-top” beam, 
defined by a step function at f = t&o as determined by the fundamental Gaus¬ 
sian in (a). Notice once again the crossing of local maxima in (c), not observed 
for unmixed Laguerre-Gaussian modes. It is this crossing that has allowed the 
beam to deliver more energy to the center of the target. 
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E. Conclusion 


Resolving arbitrary optical modes into higher-order mode spectra provides a useful tool 
for the analysis of FEL simulations and a self-consistent index of beam quality. Future work 
might include retaining the azimuthal index m in 


S{r,9,T) 


c Wo 
W 






where 




~2 

T T T 

- - — (2p -|- m -|- 1) arctan - h m9, 

Zo + Zo 


during the analysis—such a basis set would be capable of generating any beam, even those 
without radial symmetry. Free electron laser simulations have also demonstrated the poten¬ 
tial usefulness of the Hermite-Gaussian modes as an appropriate basis set when studying 
off-axis electron beams and/or resonator cavity mirror-tilt. 
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Appendix; The Parabolic Wave 
Equation 


We present the parabolic wave equation in both dimensionless and dimensional form. 
Cartesian coordinates. 

The dimensionless parabolic wave equation is 


where 


dS * ^ 

Tr = 4 ^-^' 


z 

" = V 


Vi = TTTV + TTVTT 


dx"^ dy"^ ’ 

for a propagation distance L. Using the dimensionless coordinates 




X = 



LA 


X, 


we can transform to the dimensional paraxial wave equation as follows: 


and hnally the familiar 
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dz 4 
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d^2 Qy2 

dS _ iZX fd^s d^S 

dz 471 dy"^ 


dS 


d^£ d^£ 


dz 2k V dx'^ dy'^ 


d^£ d^£ ^.d£ ^ 


Cylindrical coordinates. 
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In cylindrical coordinates we use 


-L r df \ df J 06“^ ’ 

where r = r/wi and 6 = 6, and we have 

i Id /dS\ 1 (9a^ 

4 r dr \ dr J d6‘^ ’ 
iZXll d / dS\ 1 dS^' 

dTT r dr \ dr ) r'^ d6‘^ ’ 

i [1 d r dS\ 1 dS^'l 
2k r dr \ dr ) ~^ r'^ d6‘^ ’ 

and hnally 

ld£ f d£\ 1 ( 9^2 Qg 
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